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Step meandering due to a deterministic morphological instability on vicinal surfaces during 
growth is studied. We investigate nonlinear dynamics of a step model with asymmetric step kinet- 
ics, terrace and line diffusion, by means of a multiscale analysis. We give the detailed derivation of 
the highly nonlinear evolution equation on which a brief account has been givena. Decomposing the 
model into driving and relaxational contributions, we give a profound explanation to the origin of 
the unusual divergent scaling of step meander £ ~ 1/F 1 ^ 2 (where F is the incoming atom flux). A 
careful numerical analysis indicates that a cellular structure arises where plateaus form, as opposed 
to spike-like structures reported erroneously in Ref.u. As a robust feature, the amplitude of these 
cells scales as t 1 ^ 2 , regardless of the strength of the Ehrlich-Schwoebel effect, or the presence of 
line diffusion. A simple ansatz allows to describe analytically the asymptotic regime quantitatively. 
We show also how sub-dominant terms from multiscale analysis account for the loss of up-down 
symmetry of the cellular structure. 



PACS numbers: 05.70.Ln, 68.35.Fx, 81.15.Aa 



I. INTRODUCTION 



The production of solids by Molecular Beam Epitaxy (MBE) having a surface which is abrupt on the atomic scale is 
often hampered either by a stochastic roughness or due to the presence of morphological instabilities. The stochastic 
roughness is often attributed to shot noise from the incoming deposition flux. As for determinitic instabilities, there 
are three general types of surface instabilities leading to kinetic roughnesses: step-bunching, step meandering, and 
islanding (see Fig.[j]). The two first categories are met on vicinal surfaces, while the last one can either be present on 
high symmetry surfaces, a typical mechaiusm being the Ehrlich-Schwoebel effectta, or even on a vicinal surface as a 
secondary instability of the step meanderEI. 
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FIG. 1. Schematic view of a crystal surface suffering various instabilities. Step-bunching, step-meandering, and island (or 
advacancy) formation are depicted. 
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Kinetic roughening has long been as a mystery. With regard to MBE growth on a high symmetry surface, a 
prominent example is the Kardar-Parisi-Zhangu equation introduced in an attempt to describe surface noise- induced- 
roughening, its one dimensional version is: 

d t y = aCLF + d xx y + {d x yf + 77, (1) 

where y is the surface height, and x the coordinate along the surface(|l|). Derivatives are subscripted, that is 
dty = dy/dt and so on. F is the incoming flux, a is the atomic height and fl is the atomic area. We have set 
the coefficients to unity, since only the form of the equation matters in this discussion. This equation has given 
rise to a variety of investigations both analytically and numerically (this part has included analytical treatment of 
the partial differential equation together with numerical Monte-Carlo simulations which mimics KPZ dynamics). 
Equation (|l|) is phenomenological in the sense that it is derived on the basis of symmetries. 

The KPZ nonlinearity is a natural candidate in the long wavelength limit, if dcsorption is present, or if allowance 
is made for defects (such as vacancies-often named overhangs- in the growing solid) . No derivation of that equation 
has been given so far, however. The reason is, in our opinion, the lack of a continuum description of island nucleation. 
In the absence of both defects and desorptkp (which are two usual requirements for production of solids of interest!), 
the nonlinear KPZ term is not permissiblecl. In that case the equation must have a form of a conservation law, that 
is: 

d t y = atiF - V.J, (2) 

00 that upon averaging, the mean velocity is simply be given by a£lF, as it should be; because the KPZ nonlinearity 
can not be written as a flux (i.e. as a divergence of a current) it gives an additional contribution to the grwoth 
velocity by a quantity < {d x y) 2 >^ (the symbol < .. > stands for to the average), which obviously makes no sense. 
There has thus been a variety of attempts with the aim of deriving the appropriate surface evolution equation in 
that limit. Here again no derivation from first principles is available. 

As said above, in addition to surface roughness caused by shot noise, nominal high symmetry as well well as 
vicinal surfaces, may become inherently unstableo when brought away from equilibrium. Nominal surfaces may 
develop mounds due to the ES effect. However, derivation of the appropriate surface evolution equation in that case 
is still a matter of debate, though a significant progress has been achieved. 

In contrast to nominal surfaces, vicinal surfaces in the step flwft regime have allowed to derive evolution equations 
from first principles. In a series of papers, we have showroH'El that vicinal surfaces offer a relatively tractable 
situation, though often nontrivial, where evolution equations can be extracted from basic transport and kinetic laws. 
The strategy is to first focus on derivation of step evolution equations. Once this task is achieved, it becomes then 
possible to derive the surface evolution equation. In its general form, the evolution equation is nonlocal and highly 
nonlinear. A rather simple information is extracted, however, if we focus on the long-wavelength limit: that is 
we assume that the wavelength of the step meander, and thus surface modulation, is small in comparison to the 
natural physical length (diffusion length if desorption is important, otherwise the interstep distance, which is the 
most frequent situation). More precisely the full growth equations, which are highly nonlinear and nonlocal, can be 
reduced to nonlinear partial differential equations, which are more tractable and often allow a signficant analytical 
progress as wil be shown here. n 

As shown by Bales and ZangwillE, a straight step during MBE growth may become morphologically unstable in 
the presence of an attachment asymmetry (the Ehrlich-Schwoebel (ES) effect) at the-,step. Close to the instability 
threshold, starting from the Burton-Cabrera-Frank (BCF)B model, we have shownO that the step profile in the 
presence of desorption obeys the Kuramoto-Sivahinsky equation (written in a canonical form): 

dtC = -d xx ( - d xxxx ( + (d x () 2 , (3) 

where x is the coordinate along the step (Fig.|l|), and £ designates the step position. In a similar fashion we have 
shown later that steps on a vicinal surface obey a set of coupled anisotropic Kuramoto-Sivahinsky equations!! The 
ultimate stage of surface dynamics is found to be spatiotemporal chaos. Two remarks are in order: (i) the KPZ 
nonlinearity is of the KS type -due to desorption- (ii) the first term in the KS equation has a negative sign, signalling 
an instability; there is a necessity for taking higher order derivatives into account in order to prevent arbitrary short 
wavelength modes to develop. _ 

A question of major importance arose recentlyu: if desorption is negligible, what kind of nonlinearity should one 
expect? because of the conserved character of dynamics, only terms which can be written as derivatives of a current 
are allowed. We could naively have thought that a natural candidate would be the conserved KS equation, namely 
dtC = —d xx £ — d xxxx <^ + d xx [(d x () 2 ]. A close inspection of the BCF equations, as shown here in details, reveals that 
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this is not the case, .though symmetry and conservation would dictate that form as the first plausible candidate. We 
have recently showncl that, for an in-phase train of steps, each step obeys the following nontrivial evolution equation: 



d t ( = -d x 



1 + {d x Q< 



dxC + d x 



d xx C 



(1 + (dxC) 2 ) 3/2 



(4) 



This highly nonlinear equation could not be inferred from scaling and symmetry arguments. It is related to a 
singular behavior of the amplitude of the meander that behaves as 1/F 1 / 2 when F is small. Instead of chaos, a 
regular pattern is revealed, the modulation wavelength is fixed at the very initial stages while the amplitude of the 
step deformation follows a scaling law w ~ t 1 / 2 . 

The objective of this paper is many fold. We first give an extensive derivation of the above evolution equation 
starting from the BCF model. We shall also present a general argument on why that singular behavior is present in 
the absence of desorption. A second line of the investigation concerns higher order contributions. It is clear that the 
above equation enjoys the up-down symmetry, £ — > — £. We show here that the effect of higher order contributions 
is to destroy this up-down symmetry. An important fact to be presented here is -that the step profile exhibits a 
plateau-like morphology. This contradicts the preliminary simulation given in Ref.EJ. That simulation had suffered 
from numerical inaccuracy causing spurious spikes to develop. Finally we show that though the full equation is highly 
nonlinear it has been possible to provide a quantitative analytical treatment for the step morphology, and evaluate 
the plateau width along with the meander amplitude. The results are found to be in good agreement with numerical 
results. 

This paper is organized as follow. We write down the basic equations in section O. In section [II a linear 
stability analysis is performed, which allows to evaluate the most unstable wavelength and the typical time for the 
appearance of the instability. In section IV we shall provide a general argument on the extraction of the scaling of the 
step position with the incoming flux. Section [v] presents the detailed derivation of the principal evolution equation 
(0) in the one-sided limit. We shall then present the situation where there is a finite ES barrier. Section VI deals 



with the higher order terms and their impact on the up-down symmetry. In section VII we generalize the d erivation 
of the step evolution equations to the two-sided case. Discussion and outlook are presented in section VIII. 



II. BASIC EQUATIONS 

We present the model based on that of BGEL supplemneted with asymmetric attachment kinetics as introduced by 
SchwoebelB, and line diffusion following Ref.tJ. A vicinal surface, whose mean interstep distance is £, is considered. 
On the terraces, the adatom concentration c m between steps m and m + 1 evolves according to: 

d t c m =DV 2 c m + F, (5) 

where D is the adatom diffusion constant, F is an incoming flux of adatoms from a beam, and dt denotes the time 
derivative. Once an atom is attached to the surface, it cannot detach from it (no desorption). We consider the widely 
used quasistatic limit where the concentration reaches a steady state regime on time scales much faster than that of 
step motion. We then have to solve Eq.(||) with the l.h.s. equal to zero. For implications due to non-quasisteady 
effects see Ref.EJ. 

The excursion of the m th step about its straight configuration is denoted £ m (x,t), so that its position is ml + 
£ m (x,t) + Vt, where V is the mean step velocity. We consider the case where no step overhang is present, so that 
the function is univocal. On both sides (+ and — designate the lower terrace and the upper one respectively) 
of step m, the normal diffusion flux is linearly related to departure from equilibrium with kinetic coefficients v±: 

Dd n c m \+ = v+{c m - c eq )\ + 
Dd n c m -i\- = -i/_(c m _i - c eg )\-, (6) 

where c eq is the local equilibrium concentration, and d n denotes the derivative in the direction which is normal to 
the step. More precisely d n = n.V where n = (— £ x , 1)/ \/l + is the unit vector normal to the step, and V is the 
two dimensional gradient operator: V = {d x ,d z ) where x is the coordinate along an originally straight step, and z 
the one orthogonal to it. The attachment lengths on both sides of the steps will be used later: d + = D/v + and 
d = D/v-. If c° eq is the adatom concentration close to a straight step, the concentration for a curved step is given 
by! 

c eg =c° eq (l+TK m ), (7) 
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where T = Q^/ksT (the definition of T is slightly different from that of Ref.El) with 7 the step stiffness, and n m , 
the step curvature is given by: 



[i + (^C™) 2 ] 3/2 ' 



Here for simplicity we disregard step-step elastic interaction. We shall come back to this point in the discussion. 

At the steps, mass conservation, in the limit where the adatom concentration is much smaller than that of the 
solid I/O, imposes: 

V n = ft {Dd n c m \+ - Dd nCm - X \-) + ad s [D L d s (T Km )], (9) 

where a is an atomic distance. Using Einstein's relation, the macroscopic diffusion constant along steps is defined 
as Dl = D st ac s ti where D st and c s t are the diffusion constant and the-concentration of mobile atoms along the 
step, respectively. This expression is in agreement with that of MullinstLl One may object that c s t is in factnat 
well defined along a step. We shall therefore use a more general expression derived from the Kubo formulaliiltJ: 
Dl = a 2 /tl where tl is the characteristic time for detachment of an atom from a kink. Non-equilibrium effects 
related to line diffusion are not considered in this expression. 

Two sources of nonlinearities can be identified. The first one is apparent in the boundary conditions ([?], ^|) because 
both the normal to the step and the curvature (see Eq.(||)) are nonlinear functions of the step profile. The second one 
originates the free boundary character (Stephan problem) is a hidden source of nonlinearity: the concentration field 
on a terrace -which is a nonlinear function of the position- depends on the step profile, leading thus to a nonlinear 
concentration field as a function of the step position. 

In addition to clastic interactions (not included here), steps are coupled via adatom diffusion. Dynamics are 
nonlocal in space and time. With the help of an integral formulation of the model equations, we have made explicit 
this nonlocality in a previous workEI. The use of the quasistatic approximation suppresses delay effects, whereas 
spatial nonlocality persists. 

III. LINEAR STABILITY ANALYSIS 

The linear stability analysis is the first step in any stability problem. Moreover it will allow us to prepare some 
preliminaries for the nonlinear analysis. Let us define the Fourier transform of the meander as: 

E / / CmOM) e-*"-**^*" dx dt, (10) 

where ilu is the pulsation of the perturbation of wavevector k and phase shift between two neighboring steps, <f>. 
The phase varies between and 2ir. Let us quote two special cases. The in-phase mode <f> — 0, corresponds to the 
case where all step meanders are identical, i.e. C, m {x,t) — ( m i(x,t) for any m, m! . The out of phase mode 4> = tt 
corresponds to the situation £ m (x,t) = — £ m +i(x,t). I— . 

The derivation of the full dispersion relation can be performed in this case along the same lines as in Ref.EZl. We 
shall not repeat here the calculation, but give directly the result. The quantity iuj is complex, and let us discuss 
separately the real and imaginary parts. The real part of iu> takes the form 



V \£ + d-+d + 

iy r 



(d- + d+) {q£ sinh(^) - cosh(g^) + cos(0)) + -q£ smh(q£) 



L> s J: (2 (cosh(^) - cos(0)) + q(d + + gL) smh(q£)) + aD L q 2 ] , (11) 



with q = \k\, and 



V 



V={d++ d-)q cosh(^) + (d + d-q 2 + 1) smh(q£). (12) 



Both macroscopic diffusion constants (adatom tracer diffusion constant D times coverage of mobile atoms) on the 
terraces D$ — Dilc® q and Dl along the steps enter this relation. The "bare" (tracer) diffusion constant of adatoms 
on terraces does not appear alone. 

A positive 5Re(iw) is a signature of an instability. The straight step is unstable during growth provided that a 
normal ES effect is present (d- > d+). Moreover, the most unstable mode is the in-phase mode <f> = 0. This remark 
will be exploited later. 
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The imaginary part of ioj describes propagative effects: 

Smf>>) = QFsin(4>)^(£ + d+ + (13) 

The origin of this term is quite transparent. In the limit of a straight step (q = 0), we have Qmiico) — QF sin((/>), so 
that the perturbed solution takes the form (ignoring the real part ofiw), Cm ~ e in<t>+itaF sin ^ = e i<t>(n+t(V /e) 8 in(<t>)/<l>]_ 
Here we have introduced the step velocity of the uniform train, Vb = flF£. This means that in order to travel a 
distance n£, it takes for a perturbation a time given by (n£/Vo)(<p/ sin(^>)). Since 0/sin(</>) > 1, that time is always 
longer than that needed for a uniform train (4> = 0) to travel the same distance. In other words, all perturbations 
(except the in phase one) travel forward slower that the train velocity Vq . This means that perturbations are advected 
backwards in the reference frame with velocity Vq. 



IV. SCALING ANALYSIS 



Once the instability threshold is reached, any perturbation will amplify exponentially in the course of time, so 
that nonlinear effects can no longer be disregarded. As discussed in section the set of growth equations is highly 
nonlinear and nonlocal, so that only a "brute force" numerical analysis would give a general answer. Our idea is to 
inspect the original equations and try to reduce legitimately the complexity. The key ingredient in our analysis is 
the identification of a small parameter. 



A. Scaling of space and time variables 



When inspecting the dispersion relation ( |11|) for an in phase train, one realizes that the band of unstable wavenum- 
bers extends from q = (actually this result is traced back to translational invariance; it corresponds to a global 
motion of the train) to a critical finite value q c (to be defined below). We shall assume that ql remains small in 
comparison to one, and we come back in the discussion below to the validity of this assumption. In that case Eq.(p"l|) 
takes a simpler form 

Re[iu>(q « 1,4, = 0)] = e ^ d ~^ + d J 2 (Dst + D L a)Tq\ (14) 

We have set here $ = 0, which is the exploitation of the fact that the in-phase mode is the most dangerous one. 
We consider later the situation where small deviations from the in-phase mode are taken into account. It is seen 
that the range of wavenumbers with positive iu> is given by 



f_nFff ± _^ 1/2 

\T(D s £ + D L a) 



(15) 



where f s = (d_ — d+)/(£ + d + + rf_) is a parameter describing the Ehrlich-Schwoebel effect. The demand that the 
small wavenumber expansion make a sense is satisfied by requiring e = {q c £) 2 *C 1, and this is precisely the definition 
of our small parameter 

- ^ (16) 



T(D s £ + D L a) 



This guarantees the long wavelength regime. It is important to see from the very beginning whether this limit 
is realistic, or is it rather academic. Experimental data are available on vicinal surfaces of Cu(l, 1, 17) which have 
recently revealed a meandering instability during step-flow growthliil Their data entering the expression of e which 
are best known are flF — 3 10~ 3 s~ 1 , £ — 21.7 'A. The step stiffness can be written as 7 ~ (kBT/2a)exp(Ek/kBT), 
Ek being the kink energy. From a simple "bond counting" argument, one can evaluate the adatom eouilibrium 
concentration on a vicinal surface: Qc® q ~ exp(—E a /kBT) with E a = 3Ek- Using the result of RefJia'ta from 
step fluctuations at equilibrium, we have Ek — 0.13eV. The diffusion constant on terraces takes the form D = 
a 2 UQ exp(— Ed/UbT), ^0 ~ 10 13 s _1 is an intrinsic frequency, and Ejj w 0.45eVli3. With a lattice constant of 2.55A, 
we find: D s £ = 1.4 x 10 15 | exp(-0.84eU//c B T). „ n 
Using the Kubo formulao, one can evaluate the line diffusion constant: Dl = a 2 /tl. From experimental datat3'll3: 



D L a = aD L0 cxp(-E L /k B T), with aD L0 = 6.5 x lO 18 !^ 1 and E L = 0.89eU. With these values, we find that 
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D$£/ '(Dlo) ~ 10~ 2 in the experimental temperature range. This indicates that the stabilization of steps essentially 
occurs via line diffusion in this situation. In the one-sided limit (d + = and d_ — > oo), and around 300-ftT we find 
e ~ 10~ 3 , and A c = 2ir/q c ~ I0 3 atomic distances. This result implies that a priori the longwavelenggth limit is 
appropriate. 

The active modes in the instability are those for which q£ ~ q c £ ~ e 1 / 2 , and therefore, lengthscales of interest are 
those for which x ~ IcT 1 !" 1 . The characteristic time of the instability development is given by the growth rate of the 
most unstable mode: 



2nt 



i 



T{D S £ + D L a) 



e 



-2 



(17) 



This is obtained as t m = 2n/^ie[iuj(q = q m , (j) = 0)], q m being the wavevector of the most unstable mode, related 
to q c by q m — q c /V2. This relation provides the scaling of the time variable t ~ e~ 2 . Using the above data, we find 
that the instability typically develops after a growth of a thickness of the order of 100 monolayers. 

Before proceeding further, it is instructive to analyze briefly the asynchronized train. As pointed out in an earlier 
worm, the ES effect not only induces a morphological instability of steps, but also leads to a "diffusive repulsion" 
between steps on a vicinal surface. This dynamical repulsion will force steps to evolve in-phase. The time needed 
for steps to organize in-phase in the unstable train is t$, defined as the synchronization time of the most unstable 
mode with wavevector q — q m — qc/V^ (i-e. the on e having the maximum growth rate): 

-!- ~ 9^Re[w(g,0)| g=9m ^ = o] . (18) 
t(j, 

This time corresponds to the decay of a perturbation having phase shifts of order one. From linear dispersion relation 
(Eq.p|): 



U . ( (£ + d + + d-)(Dse + D L a) 

t m ~ 6 \D s £(d- +d++ £/2) + D L a(d- + d+) 



(19) 



where t m is the typical time for the instability to develop (Eq.(jl7j)). In the one-sided limit, t^jt m ~ e. Thus, we 
expect the synchronization time to be shorter than the instability time, i.e. t$jt m <C 1 since e<l. This means 
that steps will be synchronized in the early stage of the instability. This justifies the fact to focus on small phases, 
or possibly a vanishing phase, as we shall assume later. 
For small but finite (j) we have 

m r- / ii i\i d- — d+ (. d + + d- 

Re[iu>(q « 1,0 « 1)] = - —-——±— (q£) 2 



2 i + d++d- V £ + d+ + d- 

.2 D S ,2 \ t-„2 



(Ds£ + D L a^ + IT ^-^)T^ (20) 

It is seen from the first term that for the phase shift to be relevant we must have <f> ~ ql ~ e 1 / 2 . This implies that 
the conjugate variable m (the step position along the vicinality) has the following scaling m ~ e -1 / 2 (meaning that 
one needs to travel a distance of that order to detect phase modulations). In summary we have the following scaling 
in Fourier space 

g-e 1 / 2 , cj~e 2 , <t> ~ e 1 / 2 , (21) 

and their corresponding conjugate variables in real space 

x-e- 1 / 2 , t~e- 2 , m-e- 1 / 2 . (22) 

Beside the instability character, the problem involves progative effects which are related to the imaginary part of 
iu>. Inspection of the imaginary part of the dispersion relation (|l^) in the long wavelength and small (f> limit, shows 
that Qm(iu) ~ e 3 ' 2 . This defines a fast timescale r ~ g -3 / 2 related to propagative effects -we mean faster than 
the time scale associated with the instability ~ e -2 . Since we shall mainly be interested by a synchronized train 
(in which case the imaginary part vanishes), we shall leave out this additional complication for the moment, and 
postpone this question to a forthcoming work. 
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B. Scaling of the meander 



In order to determine the nonlinear evolution equation, following our previous workB in the presence of desorption, 
we could expand all physical quantities (concentration, step position) in power series of the small parameter, the 
leading contribution would be of order e°, followed (in principle, and in a regular expansion) by e 1 / 2 , since this is 
the smallest power encountered above. This strategy worked out when desorption is present, but not in the present 
context. The reason is that the ficst contribution in the step profile turned out to be C ~ e 1 ^ 2 - This was viewed 
as an ansatz in our previous workEJ. In this paper we provide an explanation of that fact on the basis of general 
considerations. Without having resort to an explicit derivation we shall show why this scaling is inherently linked 
with the non-desorption case. 

For that purpose it is useful to identify two 'classes' of adatoms (of course just in terms of a picture): Thermal 
adatoms of concentration ct detach from a step, diffuse on terraces and re-attach to a step. Mass transport associated 
to their motion induces relaxation towards equilibrium. Freshly landed adatoms of concentration cf have not yet 
been incorporated into a step. Their attachment result in the non-equilibrium driving of the steps. We can thus split 
the full set of equations (|5j-|9j) into two pieces by writing the model equations in the following equivalent form: 

DV 2 c T = 0, (23) 
DV 2 c F = -F. (24) 

These fields obey the following boundary conditions at the steps: 

Dd n c T = ±v±{c T - c eq ), (25) 
Dd n CF — ±v±cf 1 (26) 

where the index + and — refer to both sides of the steps. They are coupled only through mass conservation at the 
steps 

V n =V F +VT, (27) 

where the driving contribution vf is proportional to the incoming flux F: 

v F = D[d n c F+ - d n c F _]. (28) 



Indeed from the equations obeyed by cf (Eqs.(|24j), (|26|)) by making the transformation cf — * cp/F one sees that F 
scales out from the equations, implying thus that cf must directly be proportional to F. 

We can extract from cp the contribution of the uniform train, which leads to a velocity given by QF£, plus another 
contribution due to step modulations which must be compatible with conservation, vf is thus the sum of the mean 
step velocity and the divergence of a flux j that describes how mass is unequally distributed between different steps, 
and different parts of each step: 

v F = {lF(£-Vj), (29) 

with j, according to what is stated above, independent of F. 

The relaxational contribution vt is a thermal part and is obviously independent of F: 

v T = [Dd n c T ]t + ad s [D L d s TK\ . (30) 

Gradients of chemical potential /i are the driving force of the relaxational contribution. Without loss of generality, 
and as long as we deal with smooth and large scale perturbations, the thermal part of the normal velocity can be 
written with help of the Cahn-HilliardEJ equation: 

ifr = V[MV4 (31) 

where M is the macroscopic mobility of the surface, and fi = Q^/n is the chemical potential. The step index m is 
omitted in this section to simplify notations. Thus we shall from now on use the scalar mobility M along x. The 
chemical potential is expressed as [i — QST/S^, where T is the step free energy. Thus, if f(d x C) is the free energy 
density, we have: 

M=-^[/W)]- (32) 
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The evolution equation of the step meander (i.e. when the step mean velocity is subtracted) now reads: 



d t C = -d x [QFj + Md xx f] . (33) 



Recall that F is proportional to e (Eq.(|16|)), so that we can set F — eF, where F is of order one. On the other 
hand x = Xe^ 1 / 2 (Eq.(|22])), so that d xx = edxx- j, M or /' only depend on derivatives of £, due to translational 
invariance. We write their argument symbolically as {d x } = {^^dx} (which is taken to mean any derivative and 



any power). Equation (33) can be rewritten as 



d t C = -e 3 / 2 d x [nFj{e^ 2 d x } + M{e 1 ' 2 d x }d X x f {e 1/2 d x }\ ■ (34) 

The central point lies in the fact that the small parameter e appears as a common factor, in the first term it stems 
from F while in the second one it originates from the second spatial derivative. 

In a regular expansion, close to the instability point we expect that the amplitude of modulation is vanishingly 
small when e«l. In reality, and this is the heart of the proof, due to the structure of the above equation, it will 
follow that no nonlinear term can enter the evolution equation, even if the amplitude were allowed to be of order 
one. There is even a stronger statement. Indeed even if C — e * H [H is of order one), with d > —1/2, we show below 
that any nonlinear term has a vanishing contribution. For that purpose we expand any function noted h (which 
represents ]■■■■) in a Taylor series 

h = h + h^'^idxH) + h 2 e 1+ ™(d x H) 2 + h.o.t., (35) 

where we have kept the smallest linear and nonlinear terms. For example a term like e® +1 dxxH <§C e° +1 / 2 dxH . 
Although our conclusion can be made at this stage, let us be more explicit. Setting T = e 2 t, Eq.(^3|) now reads: 

e a+2 d T H = -e» +2 d x [jidxH + e» +1 / 2 j 2 (d x H) 2 
+ (m + ^^khdxH) dxx (f[d x H + e §+1 ' 2 f^dxH) 2 )] + h.o.t., (36) 



Since $ > —1/2, we have e' ?+1 / 2 — > as e — * 0. Therefore, nonlinear terms are irrelevant in Eq.(p6|), and the full 
equation reduces to a linear evolution equation: 

d T H = - (jtdxxH + M f[d X xxxH) . (37) 

For nonlinearities to be relevant in Eq. ( p6 |) , we need e^+V 2 ~ 0(1), which is obtained when d = —1/2. But 
then, the expansion performed in Eq.(|35|) is a priori not legitimate. Indeed, higher order terms become relevant: 
(d x C) n ~ e"(' ? + 1 /2) _ 0(1) w hen i? = —1/2 for any integer n. We therefore expect a highly nonlinear evolution 
equation, as will be shown explicitly in the next section. 

How concentration scales with e can also be found using the decomposition of the concentration. From equations 
(|J) and @, we have c F ~ F ~ e. From Eq.(§|), (f|), and @), we find that 

OT-cl^c^-c^-c^TK-e 1 ' 2 . (38) 



Thus, u — Qcf + ct — c° eq ~ e 1 / 2 , and the concentration will be written in the following form: 

u{x,t) = e 1/2 U(x,t), (39) 

with U(x, t) ~ 0(1). Similarly, the meander will be written as: 

C(x,t) = e- x l 2 H{x,t), (40) 

where H(x,t) ~ O(l). n 

It is important to show why in the presence of desorption the expansion is regular, leading to the KS equation (El). 
With desorption the evolution equation can no longer be written in the form of a conservation law (^). Nevertheless, 
the above mentioned decomposition still holds in a slightly different form: instead of being proportional to F, the 
driving part is proportional to F — F eq , where F eq is the incoming flux at equilibrium that counterbalances ambiant 
desorption (F eq — c° eq /T, t being the characteristic residence time before desorption on a terrace). Hence, instead of 
Eq.(p3|), we have: 
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dt( = (F eq ~F)g + N d xx f, 



(41) 



where g and N are functions of the derivatives of £. The second term of the r.h.s. is now the one expected for 
non-copserved relaxation to equilibrium: it is directly proportional to the chemical potential variations (See model A 
in Ref.Ej, orE3). Linearizing this equation, we have to take g rs gd xx C since the first linear term proportional to d x C 
is a propagative term, not contributing to stabilization or destabilization (moreover, this term, not invariant under 
the x — > —x symmetry, is not allowed). We then find: 

d t ( = {(F eq - F)g + N f{] d xx Q. (42) 

The prefactor of d xx ( is the effective stiffness of the step. An instability is signaled by a negative sign of that prefactor. 
This happens when F > F c = F eq + Nof[/g. The small parameter (that is the distance from the instability threshold) 
is now e' ~ F — F c . Moreover it was found in RefJ§E3 that x ~ e' -1 / 2 and t ~ e'~ 2 . Defining as in the conserved case 
X = e n / 2 x and T = e' 2 t, and C = z' e> H, with X,T,H ~ 0(1), Eq.@ is now expanded for F w F c as: 

e'»'+ 2 d T H = -e"'+ 2 g d xx H + e' 2 ^ 1 (F F eq ) (d x H) 2 g 2 

+e m '+ 2 [N f 2 d XX {dxHf + N, f{ d x Hd xxx H] + h.o.t.. (43) 

As before, we use $' > —1/2 so that expansion ( ^5| ) makes a sense. It is seen from this equation that the leading 
nonlinear term is (F — F eq ) e' 2§ +1 (d x H) 2 . It counterbalancp, the linear term e" 5 +2 g d xx JL provided that i?' = 1. 
The nonlinear term here, is that of the Kardar-Parisi-ZhangE 2 ] and Kuramoto-Sivashinsky&El equations (see Eq.(0) 
and (||)). It could not be present in the conserved case because it cannot be written as the divergence of a flux. 
Moreover, it is non- variational, and thus it must vanish at equilibrium, as can bee seen from its prefactor F — F eq . 

It must be emphasized that the decomposition into an equilibrium an a nonequilibrium part holds in the present 
problem, but is not a general property. This does not have to be the case out of equilibrium in general (an example 
is that of step flow or sublimation in the presence of electromigration, such a decomposition between relaxation and 
driving parts has not been made possible). 



V. ONE-SIDED SYNCHRONIZED STEPS 



A. Multiscale analysis 

In addition to synchronization, we first assume for simplicity a one-sided limit (steps advance only thanks to atoms 
from the terrace which is ahead), formally defined as d + = 0, and d_ — > +oo. In this limit Eq.(^) reduces to c + = c eq 
(which is the Gibbs-Thomson condition) and dc-jdn = (atoms do not descend the steps). 

As we have shown in the last section the meander £ ~ e -1 ^ 2 , while the concentration field u ~ e 1 ^ 2 , we find it 
convenient to set £ = e _1 / 2 iJ and u = e 1 / 2 !/ , with H and U being quantities of order one. Under the assumption 
that these quantities are analytic functions of e 1 ^ 2 , we seek solutions of the form 

U = U<® + e 1 / 2 C/ (i/ 2 ) + eU W + e 3 / 2 U^ + • • • , (44) 

H = H<-V + e l/ 2 tf(l/ 2 ) + e ij(D + £3/2^(3/2) + . . . _ (45) 

In order to make explicit the e dependence, and to deal with quantities of order 1, and according to (^2|), we set: 

x = e~ 1/2 X, t = e- 2 T (46) 

It is convenient to rescale space by i and time by D/£ 2 . Performing the variable change: Z = z — ( m (x,t), mass 
conservation (j^) on terraces reads: 

= p 2 d zz U + e 1/2 (77 - 2d x H 8 XZ U - 8 XX H d z U) + e d xx U, (47) 

where p = [1 + {d x H) 2 ] 1 / 2 , an 77 = e D/QFi 2 . At the steps, the Gibbs-Thomson relation at Z = 0, and a zero-flux 
condition at Z = 1, takes the form: 

U\ z=0 = -K, (48) 
p 2 d 2 U\ z=1 = e 1 ' 2 d x H d x U \ z =i , (49) 
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where JC = fic° 9 T d X xH/p 3 . 

Mass conservation at the step (Eq.(^|)) yields 



V + e 3 / 2 d T H = p 2 d z U\ z=Q - e 1 ' 2 d x H d x U\ z=0 - e 3 / 2 8 X ( ^d x JC ) , (50) 



where /3 = D^ajDsi. The strategy is now to solve equations (^7| - |50| ) in successively higher orders in e. 

order 

To this order, Eq.(|^) reads: 

dzzl/W = 0, (51) 

which is solved by U(°) — A^Z + B^ . Equations ( |48| ) and (^) provide two conditions from which we get = 
and J3(°) = —JC^°\ No contribution to step velocity is found to th order. That is to say this order corresponds to 
the equilibrium case. 

order 1/2 

From ( fl7| ) we find that [/(V 2 ) obeys an inhomogeneous equation on terraces: 

p m dzzU <U*> = - v , (52) 

whose general solution takes the form: 

U^ = ^ V + A^Z + B<V2) . (53) 



From boundary conditions at the steps ( |48| ) and (|49| ) 



C/(V2) = B d/2) = _/c(i/2) (54) 

p (l/2)2 fe[/ (l/2)| z = 1 = dxH (0) dxB {im. (55) 

Integration constants are found to be: 

AW* = fa - a x ^ (0) M (0) ) / p (0)2 , (56) 

B (i/2) = _ iC (i/2). (57 ) 

Mass conservation at the step (|o|) determines the mean step velocity. Going back to physical variables, we find the 
expected result: V = flF£. 

order 1 

To this order, J/W obeys: 

1 



d zz U^ = [d xx H^ d z U<W + 2d x H^ d 2 xz U^ 

-2d x H^ d x H (1 ^ d zz u (1 ^ - dxxU<® 
= a + bZ, (58) 

whose general solution takes the form: 

[/(l) = ^3 + « Z 2 + A (i) 2 + B (i) (59 ) 

Once again, integration factors A^ and B^ are found from boundary conditions ( |4§| ) and (^9|): 
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\z=o = 
d z U^\ z=1 = ^ + 



d x H^ a x f/ (1/2) 

+ d x H^ dxU i°) 

- 2 d x H^ d x H^ dzUW 



(60) 



1 

1°? 



\z=i- 



Finally, mass conservation (^0|) leads to the sought after evolution equation for H^: 
M (0) = d z U^p m - d x U^d x H^ -d x ( -4m (0) 



,(o) 



(61) 



(62) 



Upon substitution of the expressions of 17 W and C/' 1 / 2 ), one realizes that terms containing /ft 1 / 2 ) cancel exactly in 
this expression, leading to a closed form for the evolution equation for H (°) : 



<9t# (0) = -d 



'x 



d x JC^ 



Going back to physical variables, we obtain: 

'QF£ 2 d x ( 



dtC = -d x 



2 (1 + (^C) 2 ) 



- (Dsi + D L a{l + (d x Cf) 1/2 ) 



d x (r K ) 



(63) 



(64) 



Besides the term proportional to Dl (line diffusion constant), this is the equation derived in Ref.B on which we have 
given a brief account. 

Introducing the step macroscopic mobility M., and the chemical potential fi = k,BTTn, the evolution equation can 
be rewritten in a more compact and enlightening form: 



dtC = -d x 



d x C, - Md s fi 



(65) 



where s is the arclength along the steps, and i± — £/[l + (c^C) 2 ] 1 ^ 2 is the distance between two neighboring steps 
measured along their normal (see Fig. || for geometrical definitions.). The effective step mobility reads: 



M = 



knT 



(66) 



The expected decomposition of step velocity (see section IV B) is clearly seen here. The first term on the r.h.s. of 
Eq.(|65|) is the driving part. To this term a simple geometrical meaning can be assigned (see Appendix The 
second term is the relaxation part with a mo bility depending on the local step orientation. Note that the present 
mobility M. and the one introduced in section IVB, noted M, differ by the scale factor [1 + (c^C) 2 ] 1 / 2 which relates 
the arc-length s to the Cartesian coordinate x. 



B. Numerical solution 

In RefEl, numerical solution of Eq.(|^) (without line diffusion) was performed using a simple Euler scheme, and 
it was found that: (i) A cellular structure takes place, the wavelength of which (the most unstable one) is fixed at 
the initial stage of the instability and no coarsening is seen, (ii) The amplitude grows like t 1 ' 2 . (iii) The shape of 
the cells is similar to the inverse error function, that is to say it develops a spike-like morphology, (iv) The meander 
is symmetric with respect to the transformation £ — > — £. It has been realized meanwhile that, though all these 
qualitative features were correct, the spikes are the result of a numerical deficiency in the original code. In order to 
cure this problem we have to resort to a special numerical treatment. Different successful attempts have been made 
but we shall here describe only the most robust numerical treatment. 

We use a powerful geometrical representation of the meandeiEj'Efl, in terms of the arclength s and the angle 9, 
oriented counterclockwise, between the normal and a given fixed direction (the z-axis direction) (see Fig.^J) . 8 is 
related to C via: tan(0) = — d x C, and the curvature simply reads: K — d s 9. 
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FIG. 2. Some definitions of the notations used in the text. 



Simple differential geometry0 provides us with the evolution equation for 9, as a function of tangential and normal 
velocities, Vt and v n : 



06 
St 



V t K ■ 



9Vn 

ds 



(67) 



Physics is invariant under a change of definition of the arclength s. This allows an arbit r ary , ,t i p e-dep endent re- 
parameterization of the curve. This "gauge" can be seen as an additional tangential velocity,E3Eja with no physical 
relevance. A particular choice that is convenient here is the one that keeps the relative arclength s/L constant in 
the course of time, where L is the total length of the curve. [This will ensure that the discretization points remain 
equally spaced along the curve. The tangential velocity readsE3E.il: 



v t = — Kv n ds' 



Kv n ds' . 

'0 JO 

The evolution equation of the meander (p4) allows one to write the step normal velocity: 



V n =~d s 



cos(6>) sin(0) 



(3 + cos(6>) 
(3+1 



d s n 



(68) 



(69) 



Dl) and spatial variables x, by \/2£/yfe, so that only one parameter 



where time t is rescaled by 4£ /e T(l Ds 
survives: = D^ajDgt. 

Derivatives along the arclength s are evaluated using a centered finite difference method. We use a backward 
differentiation scheme with variable step for time integration. This "solver" enjoys rather good precision and In- 
stability (that is to say it is unconditionally stable and optimally attenuates high-frequency (i.e. noise) components 
of the solution) that makes it well fitted to our specific problem. n 

Our present simulations show qualitatively similar behavior as that found in Ref.tl (see Fig|S|) . A maiar difference 
is revealed however: instead of spikes, the cellular structure exhibits a plateau in the extrema regionsEj, as shown 
on Fig.||. (Here we define a plateau as a region of finite slope, as opposed to regions where the slope diverges with 
time. Hence, to some scale, plateaus are curved, and this curvature does not tend to zero). The width of a "plateau" 
reaches a constant value after a transient regime. 
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Terrace diffusion 
= o 



Line diffusion 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

An x/k m 

FIG. 3. Numerical solution of Eq. (64). Meander over one period. Cells are symmetric and develop plateaus. Their amplitude 
increases with time. 



C. Analytical study 



We show here that the above numerical results can be accounted for using simple analytical arguments. We give 
here only the main results, whereas details are relegated into Appendix [b]. The central assumption is a decomposition 
in two types of regions, where two different ansatz are used. In the large slope regions, a multiplicative variable 
separation is used: 



Cs(x,t) = A(t)g(x), 
while and additive variable separation is performed for the plateau regions: 

( p (x,t)=B(t) + h(x). 



(70) 



(71) 



An additional constraint coming from mass conservation allows to determine quantitatively the asymptotic behavior 
by matching these two solutions. The amplitude of the meander is found to behave as: 



Cmax Cmin — 2gq^ ^ ; 



(72) 



where ao is calculated in Appendix [b|. The rescaled meander ((x,t) ft 1 / 2 converges to a well defined profile, and 
looks as if plateaus were formed in the extrema regions. The width of a plateau is defined as Aq/2. We find: 



A = \J((3), 



(73) 



where j3 = Dsl/(Ds£ + D^a), X c = 2ir/q c is the largest stable wavelength from linear analysis, and / is a function 
given in Appendix |b|. Since I((3) decreases monotonously from 1(0) = 1 to 1(1) w 0.54, the plateau size increases 
as line diffusion is increased, and is always smaller than A c /2. Good quantitative agreement between the numerical 
solution of Eq. (^34j) and these analytical predictions is found (see Figji ) . 
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0.2 0.4 0.6 0.8 1 

P 

FIG. 4. One sided case: The prediction for the normalized plateau size I(/3) (solid line) is compared with numerical solution 
of Eq.(64) (symbols), which corresponds to the one-sided limit. 



VI. FRONT-BACK SYMMETRY BREAKING 



The expansion performed in section V A can be pushed to next order following the same strategy. We shall here 
merely give the result and details can be found in Ref.ta. Instead of a closed equation for H^°\ here two coupled 
dynamical equations for if' ' and f/^ 1 / 2 ) are obtained. Going back to the physical quantity £ = e~ l ' 2 {H^ ' + 
e 1 /2jj(i/2)^ C0U pled equations can be recast into a single equation for Q: 



dtC 



d_ 

dx 



±d x ( 1 



~3~ 



2l x 



2) 



ds 



where the macroscopic mobility of the step reads: 

D s £± + D L a D s I 2 k 



2k B T 



(74) 



(75) 



Hence, to this order, correction to Eq.(|64|) are proportional to step curvature. 

As before, a geometrical formulation with rescaled time and space variables, is used. We now have two parameters 
(3 and e in the normal velocity: 



v n = -d s 



roM«) ^i^if9) + ( ft^^P j d s n-^~6K f | (cos(20) + 2) sin(0) + d s K 



(76) 



The same qualitative features as for Eq.(f34|) are observed: the wavelength is fixed at early stages by the one corre- 
sponding to the fastest growing mode, and the step roughness increases with time as t 1 / 2 (see Figj^). The interesting 
fact is the symmetry of the shape: the cells do not enjoy the up-down symmetry, ( — > — £. Clearly, Eq. (|64|) is 
invariant under the up-down symmetry £ — > — £. The symmetry breaking originates from the new terms as can be 
seen by performing the transformation £ — > — £ on Eq. (f74|) . As shown in Appendix |^, these terms affect the relative 
sizes of the back and front plateaus, but the t 1 ! 2 scaling law for the roughness amplitude seems to persist as a robust 
feature. 

The,Eesults obtained in this section arc in complete agreement with full simulations based on a solid-on-solid model 
m RefJI Hence we have succeeded in extracting the relevant dynamics of step meander by means of the multiscale 
analysis. 
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Terrace diffusion 
p = o 



Line diffusion 
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xA m xA m 

FIG. 5. Meander evolution when subdominant contributions (to order 1/2) are taken into account: front-back symmetry is 
broken. 
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FIG. 6. One sided case: The roughness of the meander obeys a scaling law: w ~ t 



1/2 



VII. TWO-SIDED STEPS IN PHASE 



In the two-sided regime (i.e. d+ and gL finite), a similar multi-scale analysis can be performed for in-phase steps. 
We find: 



d t ( = d x 



OF £j(d--d + ) 1 

2 x( 'd + + d- + e± (i + ^c) 2 ) 1 / 2 



D L a + D t 



£ 2 + £ ± (d + + d-) 
d+ + d-+£± 



d x (T K ) 



(77) 



Although this equation looks more complicated, the meander evolution is qualitatively similar to that found in the 
one-sided case (which is recovered by taking the limit d_ — > oo in Eq.(|77])). Indeed, plateau formation and power 
law behavior of the roughness (with the same exponent ~ i 1 / 2 ) are also found in the two-sided case. 
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More details on step behavior, such as the plateau size, can be gained from the analytical investigation of Eq.(77), 
as shown in Appendix |b[ In the pure line diffusion regime D^a 3> D$l, we have: 



A = XJ(S), 



(78) 



where Ao/2 is the plateau size, 5 = £/(£ + d + + and I is the same function as in Eq.(|7^). In the pure terrace 
diffusion case D^a <C Dgl, we find: 



A = A c /(1 - 5). 

These result are in good agreement with numerical solution of Eq.(|7^) (see Fig]?]). 

1.4 



I( P ): Analytical results 
Numerical results: terrace diffusion 
-Numerical results: line diffusion 




(79) 
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5 , (1-3) 

FIG. 7. In solid line is plotted the normalized plateau size 7(/3). Symbols represent numerical results: y refers to the pure 
terrace diffusion case, and * to pure line diffusion. 



VIII. DISCUSSION AND SUMMARY 



Starting from the BCF model, we have extracted a nonlinear evolution equation for the step meander. This 
equation is highly nonlinear, and thus, could not be expected from traditional phenomenological approaches, where 
linear terms are simply supplemented with an additive nonlinear term, as in the case of KS or KPZ equations. 

A central result of the present study is that the late time power law behavior of the amplitude of the meander 
w ~ t 1 / 2 is a robust feature, regardless of the details of the evolution equation (one-sided, two-sided, line diffusion....). 
It is an important task for future investigations to see whether this property could be derived directly from the basic 
BCF equations. A second interesting point which is worth of mention is that higher order terms destroy the up- 
down symmetry, but the power law Wr-r^ t 1 / 2 remains unaffected. The same conclusion follows from full lattice gas 
simulations as briefly reported in Ref.u. A natural question arises: could the amplitude temporal increase continue 
to evolve without bound in all circumstances until the surface breaks up into a lamellar-like pattern, or is there a 
physical mechanism, not accounted for here, leading to saturation of the amplitude? Experimental observation of 
this instability^ seems to show such a saturation for the case of Cu(l,l,17), while experiments 2 !! on Si(OOl) does 
not reveal a hint towards a saturation. Possible candidates for amplitude saturations are (i) strong anisotropy, (ii) 
elastic step interactions. We hope to report along these lines in the near future. 

It is worth pointing out that the use of equilibrium formula to evaluate the stabilizing line diffusion effect could be 
criticized, since densities of kinks and of mobile atoms along steps depend on growth conditionsES. With regards to 
line mobility, our analysis allows extraction of the geometry dependence of the mobility for large meander amplitude. 
This treatment should serve as a basis for the nonlinear study of relaxation towards equilibrium (i.e. thermal 
smoothening) of large perturbations on vicinal surfacescj. 

Perhaps one of the most striking result is the manifestation of rather 'stringent' plateaus, which are likely linked 
to the non-standard character of the evolution equation. It is interesting to note that the plateaus are a feature of 
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a continuum theory, a finding which is to be compared to a long standing problem in the context of ES-induced 
mound formation^. In all previous studies, mound plateaus were indeed considered as a signature of the breakdown 
of continuum theories. We have shown here, in contrast, that a single equation in the continuum limit can produce 
such plateaus, without having resort to specific ingredients in the angular region. It is not yet clear what kind of 
equation in the continuum limit would describe these dynamics for mound formation. Is it similar or not to the one 
encountered here? These questions constitute an important line for future inquiries. 
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APPENDIX A: GEOMETRICAL ORIGIN OF THE DESTABILIZING TERM 



In Eq.(p4|), the relaxation term is interpreted as a Cahn-Hilliard contribution. We present here a derivation of the 
destabilizing term from geometrical considerations in the one-sided limit. 

Let us consider a curved part of the step as shown in FigJ|. In the one-sided model, step motion results from 
incorporation of adatoms from the terrace ahead of it. Mass conservation for an element of terrace surface AS, 
hatched in FigM reads: 



vAx = QFAS + jj_ (x) - j± {x + Ax) 



(Al) 



where v is the step velocity along the z axis, and Air the extent of the step element CC along the x axis. The 
number of atoms entering the step is v n As — vAx. j±(x) is the total flux accross the BC segment in Fig.||. AS* is 
written as: 

AS pttAx- A{x)+A{x + Ax), (A2) 
where v is the step velocity along the z axis, and A{x), the area of the triangle ABC on Fig.|[ is a function of d x C,: 

I 2 d x ( 



A{x) 



■ cos((9) sin(0) 



2 1 + (cU) 2 ' 



(A3) 



where 9 is the angle between the z axis and the normal to the step. In the long wavelength limit, the local geometry 
of the terrace is completely described by £±, the length of BC on Fig.||, k, the step curvature, and their derivatives 
with repect to the arclength s along the steps. Since the flux j± can only come from a variation of the local geometry 
along s, we have, at most 



j± ~ d.l±. ~ d xx ( .4 - <), C , 



(A4) 



whic h sh ows that terms containing j±_ can be neglected to leading order in Eq.(Al) . Combining Eq.(Al), Eq.([A2|), 
and (A3), and letting Ax going to zero, we find: 



(Q.F1 2 d x ( 



V 2 l + ^C) 2 , 

Once the mean step velocity V — ClFl is subtracted, we recover the first term of Eq.(|64|). 



(A5) 




Ax 



FIG. 8. The element of terrace area AS that feeds one element of step, is hatched. Its area (approximated by that of 
BB'C'C) is the sum of £Ax, the area of AA'C'C, minus A(x) is the area of the triangle ABC, plus A(x + Ax) the area of 
A'B'C. 
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APPENDIX B: LATE TIME BEHAVIOR 



In this appendix we derive analytically the main results obtained numerically. Despite the highly nonlinear 
character of the evolution equation, we show here that some simple ansatz allows to describe the asymptotic regime 
with good accuracy. 



1. Large slope regions and extrema regions 



We found the conserved evolution equation of the step meander: 

dt< = -d x j[Q, 

with the mass flux (see Eq.(^7|)): 



(Bl) 



m = 



(i + {d x Q 2 ) 1 ' 2 



« + (l + (0 x C) 2 ) 1/2 



D L a + D s l 



l + 5(l + (c«) 2 ) 1/2 



3 



8 + (l + {d x 2 Y' 2 J xx \(i + (d x Q 2 y/ 2 



(B2) 



where S = l/(d+ + and a — HF£ 2 (d- — d+)/2(d+ + For the large slope region we make use of the variable 
separation: 



( s {x,t)=A{t)g(x), 



(B3) 



where A > 1, and d x g =/= for any value of x. Substituting in Eq.flBSp, one finds that the destabilizing term 
(proportional to a) dominates and: 



AA' = a^-r = C, 
99 



(B4) 



where C is a constant, and the prime stands for the derivative. The late time solution of these equations reads: 

A = (2Ci) 1/2 (B5) 
g = (2a/C) 1/2 err 1 (4x/\ s ) . (B6) 



A s being a constant, and erf(a;) the error function. Inserting these expressions in Eq.(B3), we find that the meander 
does not depend on C: 



Ux,t) = 2(ai) 1/2 err 1 (4a;/A s ) , 



(B7) 
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FIG. 9. Asymptotic meander morphology. A m is the period of the meander and is also the most unstable wavelength 
obtained from linear stability analysis, ho is the amplitude of the meander and Ao/2 is the width of a plateau. In the large 
slope regions (of width As/2), the meander can be fitted by an erf like function (see Eq.(B7)). jo is the mass flux coming from 
these large slope regions toward the plateaus. 



This solution describes regions of large slopes, but is not expected to accurately describe the shape around the 
extrema of C where the slope d x C approaches zero. In those regions of width Ao/2 and of meander amplitude of the 
order of ho (see Fig.ffl, global mass conservation implies: 



2 jo = y^/iq, 



where jo is the mass flux coming from the large slope regions. From Eq.(B2), we have: 

Jo 



(B8) 



(B9) 



where xq — (A m — Ao)/4 is the abscissa of the crossover point between high slope and extrema regions, and the 
peri od o f the meander A m is that of the most unstable mode obtained from linear analysis. Using ho = £ s (xo,t) and 
Eq.@: 



Ao 
4 



A,: 



-erf 



h 



so that Eq.(B8) now reads: 



1 a^A, 



= d t h 



■ exp 



A m A s 

~a T erf 

4 4 



2(crf) 1 /2 



fro 
2(crf) 1 /2 

h a 



(BIO) 



2(crf) 1 /2 



(Bll) 



This equation has the trivial solution: ho = ciot 1 / 2 . Using this solution and Eq.(BlO), Ao is seen not to depend on 
time. In the extrema region, we therefore look for solutions of the form: 



C P (x,t) = B ± (t) + h(x). 



(B12) 



with B±(t) = iaot 1 ' 2 , where the plus an d m inus signs refer to the maxima and the minima regions respectively. Upon 
substitution in the evolution equation (B2), we find that the problem amounts to finding the stationary solutions 
d x j[h(x)] = 0. Looking for solutions with left-right symmetry x — > —x we finally have to solve 



j[h{x)]=0. 

This will be exploited in the next section. 

The parameters A s , Aq and ao are not independent. Usin g Ec 



(B13) 



A s and ao can be determined as a function of Ao- From Eq.(Bll 

a 



( |B8| ) and Eq.(BlO), we have two relations, so that 
), we get an implicit equation for ao: 



A„ 

An" 



- 1 = 



2aV2 



exp 



erf 



ao 



2a 1 /2 



The expression for A s is: 



A * = A °v^(^)exp (^) 



(B14) 



(B15) 



Hence, the asymptotic behavior of the meander (defined by Eq.(B7) for high slopes and Eq.(B12) for small ones) 
only depends on the size of the extrema region Ao/2 (since the two parameters ao and X s are linked to Ao). How the 
plateau size is related to the model parameters will be considered in the next section. 
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2. Plateau size 



In the one-sided limit we have (5 = 0, and Eq.(B13) yields (in view of Eq.(B2)): 



= a 



ti 



1 + h' 2 



D L aT 



Introducing the abbreviation: 



DsW 

1 + h' 2 + (1 + h' 2 ) 1 / 2 



h' 



b! 



(1 + ^2)1/2 



(l + /l' 2 )V2 



we can rewrite it in a familiar form: 

np s e + D L a)T 

\ a 



/3(1 - m 2 ) 1 / 2 + (1 - fj) 



dU 
dm ' 



(B16) 



(B17) 



(B18) 



analogous to that describing the motion of a particle of position to as a function of time x in a potential U. We have 
defined (3 = D s i/(D s t + D L a), and: 



C/(m) 



to dm 



o f3(l - m 2 ) 1 / 2 + (1 - (3) 



(B19) 



Multiplying Eq.(B18) by to' and integrating with respect to x, we get the analogue of the "energy conservation" 
condition: 



1 f(D s £ + D L a)T 
2 



to' 2 + U{m) = U{m ), 



(B20) 



where mo is the turning point (i.e. to' = when to = too). We look for solutions having two "vertical" tangents 
where the step slope diverges: d x £ — > ±oo, (i.e. mo — * ±1) in order to match the extrema solution with the large 
slope region. The size of the extrema region reads: 



Ao 
2 



A /2 



dx 



mo dm 
- m „ ml 



where toq — > 1. Using Eq.(B20), we find: 



where 



Ao 

Am 



2tt 2 



, V2 ' 
(D s e + D L a)T^ 1/2 



is the the wavelength of the most unstable mode obtained from linear analysis, and 



UP) 



dm 



irV2 7_i [C/(l)- C/(to)]!/2 



(B21) 



(B22) 



(B23) 



(B24) 



is plotted in FigJ^. 10) is a decreasing function with 1(0) = 1 and 1(1) w 0.54. Hence the extrema region size is 
finite and always smaller than A c = X m /y2. 

The meander variation in this region, as compared to the total amplitude of the meander, decreases as t~ 1 / 2 , and 
the step looks as if plateaus were present. 

The reader is invited to repeat the calculation in the two-sided case -where S is finite, in presence of pure line 
(j3 = 0) or terrace diffusion (j3 = 1). Surprisingly, the same integral / appears. Let us define 6 = £/(£ + d + + 
We find in the pure line diffusion case: 
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Ao_ 
A m 



and in the pure terrace diffusion case: 



Am 



V2- 



V2 ' 



(B25) 



(B26) 



Hence, A never exceeds X m /y/2 = A c the largest wavelength for which the meander is linearly stable. 

We can use a similar treatment to analyze the case with higher order terms (eq.^ij) . The main point is that terms 
proportional to e 1 / 2 do not affect the long time behaviour obtained from the ansatz |B3| ). Consequently, in large 
slope regions, we expect once again £ ~ t 1 / 2 . Terms breaking the front-back symmetry will affect differently maxima 
and minima regions, as seen in Fig.^J, because the effective potential U(m) is not invariant under the m — * — m 
transformation anymore, which corresponds to the up-down z — ► — z for step meander. We shall not develop here 
further this point; for more detail s se ecj. 

The numerical solution of Eq.(Bl) is performed in order to check the validity of the analytical results. First, 
the qualitative profile of the meander is in good agreement with the above description. We found the predicted 
scaling of the amplitude of the meander ~ t 1 ^ 2 in all simulations performed so far, expect in the case of very small 
kinetic lengths (d+ +d-)/l < 10~ 2 , where we were not able to explore the late time behavior, due to bad numerical 
convergence. 

As shown in Fig. 4 and 7, the observed plateau size, is in very good agreement with the prediction of Eq.( B22] , 
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FIG. 10. The solid line represents the prediction for the amplitude ao as a function of xo/X m - The dashed curve is the 

prediction of the dimensionless ratio A s /A m . Symbols represent data from numerical solution of Eq. (64) (the one-sided case), 
as is varied. 



The value of ao is extracted from the evolution of the meander amplitude via the relation: 
A s is calculated from a fit of d x ( at ( = 0: 



d x Cs\c=o 



^{atf/ 2 . 
A., 



(B27) 



(B28) 



In Fig. 10, both nu meric al values are compared to the prediction of Eq.(B14) in the one-sided limit, where Ao is 
calculated from Eq.( |B22| ). Once again, good agreement is found. 
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